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The calculation of two- and four-particle observables is addressed within the framework of the 
truncated polynomial expansion method (TPEM). The TPEM replaces the exact diagonalization 
of the one-electron sector in models for fermions coupled to classical fields such as those used 
in manganites and diluted magnetic semiconductors. The computational cost of the algorithm 
is O(N) - with N the number of lattice sites - for the TPEM which should be contrasted with 
the computational cost of the diagonalization technique that scales as 0(N 4 ). By means of the 
TPEM, the density of states, spectral function and optical conductivity of a double-exchange model 
for manganites are calculated on large lattices and compared to previous results and experimental 
measurements. The ferromagnetic metal becomes an insulator by increasing the direct exchange 
coupling that competes with the double exchange mechanism. This metal-insulator transition is 
investigated in three dimensions. 
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I. INTRODUCTION 

Strongly correlated electron systems continue to be 
one the most important areas of research in condensed 
matter physics. In this context, spin-fermion models, i. 
e. fermion systems coupled to classical fields, provide a 
bridge between the full Hubbard model and models with 
only classical degrees of freedom. The later fail to cap- 
ture important physical aspects, like carrier mediated fer- 
romagnctism in manganites while the former are usually 
too complex to study without approximations due to the 
sign problem^ 

Spin-fermion models contain all the physical proper- 
ties of a strongly correlated electron system, including 
the effect of the carriers and inhomogeneous local charge 
density. They have been successful in the description of 
many materials, e. g. manganites^ and diluted magnetic 
semiconductors (DMS) £ Inhomogeneous phases in these 
materials have been found to be of importance in ex- 
plaining many physical properties. It was also found that 
mean field theories are insufficient to obtain these phases 
in DMS systems^ but that information from these mod- 
els must be extracted without approximations. In 2005, 
a spin-fermion model for multiferroic materials was also 
studied with the diagonalization method^ 

The traditional method to calculate observables has 
been to integrate the fermions exactly at a finite temper- 
ature for every configuration of classical fields and then 
integrate these fields by means of a Monte Carlo proce- 
dure. Although this method - which will be referred to 
as DIAG - is accurate and has lead to many advances, 
e. g. the theoretical observation of phase separation in 
manganites^ and clustered states in DMS^ it has limi- 
tations. In fact, the largest systems that can be studied 
with this procedure are lattices with less than approxi- 
mately 300 sites. 

In 1999, Furukawa and coworkers^ introduced a poly- 
nomial expansion of the fermionic density of states as a 



way of avoiding the diagonalization of the fermionic sec- 
tor. This polynomial expansion method (PEM) scales as 
0(N 3 ) as opposed to the diagonalization method that 
scales as 0(N A ), N being the size of the system. The 
method could also be parallelized, which is practically 
impossible for the diagonalization method. In 2003, the 
same researchers^ developed a truncated version of the 
PEM, that will be abbreviated by TPEM that scales lin- 
early, i. e. O(N), with the size of the system. All these 
methods do not introduce systematic errors or uncon- 
trolled approximations and the results obtained converge 
to the exact result as the expansion cut-off tends to in- 
finity. 

Although the TPEM is becoming a standard way to 
avoid the diagonalization^SiiSiiiiiS observables that in- 
volve two- and four- fermionic operators have not been 
studied systematically. Here, it will be shown how to ef- 
ficiently calculate these observables and the formalism is 
applied to a double exchange model for manganites. The 
emphasis will be on the calculation of the spectral func- 
tion, A{k,uj), and optical conductivity, cr(w), although 
the treatment applies to other observables as well. 

Another algorithm that improves the integration of the 
one-electron sector in double exchange models is the "Hy- 
brid Monte Carlo" (HMC) algorithm^ This method uses 
the path integral representation, introducing imaginary 
time and a HMC procedure to evolve the system. The 
TPEM works better than the HMC at low temperatures 
where the HMC presents increasing computational cost 
due to the time discretization. Another reason why we 
have preferred the TPEM is that the HMC method is 
applicable when the bands are connected and do not ex- 
tend over a wide range of energies, as is the case of finite 
coupling systems. The TPEM also allows for easy paral- 
lelization, improving the performance even more. 

The contents of this paper is the following. In £111 Al 
a review of the TPEM for spin fermion models as in- 
troduced originally by Furukawa et al. is presented. In 
EIII Bill CI the formalism to calculate two-particle observ- 
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ables is discussed as well as examples and in particular 
the formula for the spectral function, A(k,u>). In Ull Dl 
the formalism is extended to include four-particle obscrv- 
ables and to allow for the calculation of the optical con- 
ductivity. The parallclization of the algorithm and its 
complexity is discussed in £111 El Then, in Section IIIII 
these observables are calculated for a model of mangan- 
ites. The metal to insulator transition is studied in three 
dimensions. Finally, in Section llVI we present the most 
important conclusions and discuss further uses of these 
techniques. 

II. THEORY 
A. Review of the TPEM 

The brief review presented here follows closely Refer- 
ence Consider a model defined by a general Hamil- 
tonian, H = Y^ij a p Ha^ia j/3 (0)^/3, where the indices 
i and j denote a spatial index and a and (3 are inter- 
nal degree(s) of freedom, e. g. spin and/or orbital. The 
Hamiltonian matrix depends on the configuration of one 
or more classical fields, represented by <p. Although no 
explicit indices will be used, the field (s) 4> are supposed to 
have a spatial dependence. The partition function for this 
Hamiltonian is given by: Z = J d<p Sn( n l ex P(~/^(0) + 
(3[iN)\n) where \n) are the eigenvectors of the one- 
electron sector. Here (3 = l/(ksT) is the inverse temper- 
ature. The number of particles (operator N) is adjusted 
via the chemical potential \x. The procedure to calculate 
observables (energy, density, action, etc) is the following. 
First, an arbitrary configuration of classical fields <fi is 
selected as a starting point. The Boltzmann weight or 
action of that configuration, S(<p), is calculated by diag- 
onalizing the onc-clcctron sector at finite temperature. 
Then a small local change of the field configuration is 
proposed, so that the new configuration is denoted by <p' 
and its action is evaluated to obtain the difference in ac- 
tion AS = S((f>') — S((f>). Finally, this new configuration 
is accepted or rejected based on a Monte Carlo algorithm 
like Metropolis or heat bath and the cycle starts again. 
In summary, the observables are traditionally obtained 
using exact diagonalization of the one-electron sector for 
every classical field configuration and Monte Carlo in- 
tegration for those classical fields^. The PEM/TPEM 
replaces the diagonalization for a polynomial expansion 
and is briefly d escrib ed here but further details can be 
found in Ref. F^ l8ll5i These definitions will be helpful 
later in deriving formulas for many-body observables. 

It will be assumed that the Hamiltonian H((f>) is nor- 
malized, which simply implies a re-scaling: 

H -> (H-b)/a 

b = (E max + E min )/2, (1) 
where E max and E m i n are the maximum and minimum 



eigenvalues of the original Hamiltonian respectively. This 
in turn implies that the normalized Hamiltonian has 
eigenvalues e v e [—1,1]. 

Observables can be divided in three categories: (i) 
those that do not depend directly on fermionic opera- 
tors, e. g. the magnetization, susceptibility and classical 
spin-spin correlations in the double exchange model, (ii) 
those for which a function F(x) exists such that they can 
be written as: 

A(tj>)=j F(x)D(cf),x)dx, (2) 

where £>(</>, e) = ^ v <5(e(^) — and t v are the eigenval- 
ues of H(4>), i- e. D((j),x) is the density of states of the 
system. Finally, category (iii) would include fermionic 
observables that cannot be expressed as in Eq. (J2J), in 
particular two and four-particle observables. In this pa- 
per the emphasis will be on this last category, that will 
be analyzed in the remaining sections, since dynamical 
observables are expressed in terms of two-particle observ- 
ables, e. g. the spectral function and four-particle linear 
responses, e. g. optical conductivity, that had not been 
studied before with the TPEM or on large systems. 

For category (i), the calculation is straightforward and 
simply involves an average over Monte Carlo configura- 
tions. Category (ii) includes the effective action or gen- 
eralized Boltzmann weight and this quantity is particu- 
larly important because it is calculated very frequently 
throughout the Monte Carlo procedure to integrate the 
classical fields. Therefore, we first briefly review the cal- 
culation of these types of observables. As Furukawa et 
al. showed, it is convenient to expand a function F (x) in 
terms of Chcbyshcv polynomials in the following way: 

oo 

F(x) = f m T m (x), (3) 

where T m (x) is the m— th Chcbyshcv polynomial evalu- 
ated at x. Let a m = 2 — <5 m ,o- The coefficients f m are 
calculated with the formula: 

fm. = J a m F(x)T m {x)/(TrVl-x 2 )- (4) 

The moments of the density of states are defined by: 

Mm(<£) = Y^=r( v \ T ™{H{4>))\v), where N dlm is the di- 
mension of the one-electron sector. Then, the observable 
corresponding to the function F(x), can be calculated by: 
A{<l>) = Em fmfJ>m{4>)- I n practice, the sum in this equa- 
tion is performed up to a certain cutoff value M , without 
appreciable loss in accuracyAiii The calculation of \i m is 
carried out recursively. \v\m) = T m (H)\v) = 2H\v;m — 
1) — \v;m — 2) and hence fi^m = ^2 v {{ m \ v \ v \ m ) ~ 1) 
and fi-zm+i = J2u(( m > v \ v 'i m + 1) — (v; 0|i>; 1)), are used 
to calculate the moments. The process involves a sparse 
matrix-vector product, e. g. in T m (H)\v), which, when 
considering the trace operation, yields a cost of 0(N 2 ) 
for each configuration, implying 0(N 3 ) for each Monte 
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Carlo step. In addition, an improvement of the present 
method has been proposed^ based on a truncation of the 
matrix-vector product mentioned before and it turns out 
that the resulting algorithm has a complexity linear in N. 
This approximation is controlled by the small parameter 

The function corresponding to the effective action for 
a configuration <j> is defined by F s (x) = — log(l + 
exp(— (3{x — //))) and S((j>) admits an expansion as before 
with coefficients corresponding to F s (x). This ob- 
servable is calculated very frequently in the Monte Carlo 
integration procedure and so it is important to calculate 
it efficiently. In actuality, only the difference in action, 
AS = S ((/>') — S(<f>) has to be computed at every change 
of classical fields. The authors of Ref. have also de- 
veloped a truncation procedure for this trace operation 
controlled by a small parameter, et r and the current value 
of the effective action can be calculated faster using the 
moments of the expansion from the previous MC update. 
This truncation is based on the observation that if <f> and 
cf>' differ only in very few sites (as is the case with local 
Monte Carlo updates) then the corresponding moments 
will differ only for certain indices. More details can be 
found in Ref. la- 
in what follows, the size of the Hilbert space will be 
denoted by Ndim and it will depend on the size of the 
lattice as well as on the particular model to be studied. 
For a one-orbital double exchange model on a lattice of 
N sites and finite coupling, Ndim = 2iV; the factor of 2 
is due to the spin degree of freedom. 



B. Density-Of-States 

The fermionic density of states of the system for a 
configuration of classical fields <fr is given by N^(u>) = 
^ A 5(w' — 6a), where lj' = (u> — b)/a. Then, the corre- 
sponding function for D^(u) in Eq. (J2J is F(x) = 5(ui'—x) 
and the corresponding fixed coefficients are calculated 
from Eq. The end result is 



C. Two-particle observables 

The previous derivations can be slightly modified to al- 
low for the calculation of two-particle observables. These 
can be written as: 



A ia ,j a >{4>)= I F(x)D iada ,((j),x)dx, (6) 



where 



Di, 



x) 



E 

A 



u i, ia u 3 a '> xS ( x - eA ) 



(7) 



and F(x) is a function that defines the observable in ques- 
tion. The quantities Uj a> \ and e\ denote the eigenvectors 
and eigenvalues of the normalized one-electron Hamilto- 
nian, H (<f>), respectively. (In general, i and j are spa- 
tial indices and a and a' are indices corresponding to 
internal degrees of freedom if any. Spatial indices will 
always be denoted without arrows or bold letters inde- 
pendently of the dimension since on any lattice a vector 
v can be mapped into a single integer. Also the notation 
i + j is meant to represent the lattice site given by the 
vectorial sum of the vectors corresponding to i and j re- 
spectively.) For example, in the context of a one-orbital 
spin-fermion model with spin, consider the observable: 
(4<*0 - Ea Ul i(T U ja , x f(/3(ae x + b - M )), with f(x) = 
^ipjj and the F(x) of Eq. © is F(x) = f(fi(ax + b- fi)) 
for this particular observable. (The normalization factors 
a and b where defined in Eq. (JJJ.) Now, F(x) is expanded 
in Chebyshev polynomials, F(x) = ^™ =0 /m^mW; an d 
the coefficients f m are calculated by inverting this for- 
mula. In all equations (■■•) denotes the average over 
the fermionic sector for a fixed configuration of classical 
fields. The moments of Di a _j a i arc defined by: 

A*ia ,ja';m {4>) = (ia\T m (H(cb))\ja'), (8) 

and therefore: 

Aict ) jct / (0) = ^ fmf^i<x,ja' ;m (0) j (9) 



(5) 



A Monte Carlo average over classical fields is generally 
needed to obtain the fermionic density of states, i. e., 
N(co) = (iV^(w))^, and for that it suffices to calculate 
{Hm{4>)) 4>> as can be deduced from the form of the last 
equation. 

The sum in Eq. O) is truncated to a certain cut-off 
M (see, e. g. Ref. llRTI and such an abrupt truncation 
generally results in unwanted Gibbs oscillations. Silver 
et al. have considerecU&ii instead smooth truncations, 
where the moments are multiplied by damping factors, 
g rn . A possible choice of the damping factors, g m , can be 
found in Ref. ITIT7I 



which is the formula used for the calculation of the 
observable Ai a j a > with the TPEM. The computational 
complexity of the calculation of iMaja'ym can be inferred 
from the computational complexity of the calculation of 
the moments of the density of states^ and is O(l) for 
TPEM and 0{N) for PEM. In general, non-local observ- 
ables involve a sum over the index i or j in the previ- 
ous equation. For instance, the kinetic energy is writ- 
ten as K = —tJ2(ij) a( c l<j c jo- + H.c.) and, therefore, the 
complexities are multiplied by a factor of N. However, 
this still gives O(N) complexity when the truncated ex- 
pansion is used. It is important to remind the reader 
that two-particle observables are not required at every 
change of classical fields in a Monte Carlo simulation 
but are computed at most only once per Monte Carlo 
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step. One observable of particular interest is the spec- 
tral function, A(k,oS) that can measured with Angle Re- 
solved Photoemission Spectroscopy (ARPES). This ob- 
servable can be expressed in terms of the eigenvalues e\ 
and eigenvectors U\^ a of the normalized Hamiltonian^: 

A-4>{r,u)) = Y,i, a ,\ u \,ia(4>)U{i^) a ,\{<t>)5{u' - ex), where 
lu' = (u> — b)/a is the normalized frequency. Then, it is 
straightforward to show that the function F(x) of Eq. © 
for A(r,uj) is F A (x) = S(u>' — x), which allows for the 
calculation of the Chebyshcv coefficients in terms of the 



moments of D 



A(r,w) 



MT,jcrV 



with the help of Eq. Q : 



r)a\m 



irVl 



■ iv'- 



(10) 



(Note that a m — 2 — S m .o was defined previously. 



and the elements (icr\j x \ja') can be easily calculated: 
they are equal to ±ite5 a _ a i if j = i ± x or otherwise. 
Now it is possible to calculate, rjmn{4>) for any classical 
field configuration <f>. Finally, a Monte Carlo averaging 
over the configurations <j> is performed that yields u(uj): 

1 f°° 

<t(u) = - y^fann(^)) <pT m (x)T n (x + Lj) X 

U J -™ ran 

(f(/3x)-f((3(x + cj)))dx. (14) 

Note that the Monte Carlo procedure has to average only 
M x M quantities (r]mn(4>)}<l> and the computation of a{u>) 
using the previous equation needs to be performed only 
once per Monte Carlo step. It is also worth remarking 
that our method differs from the one proposed by Iitaka 
et aliSS because it does not introduce random vectors and 
it is applicable at finite temperature. 



D. Four- particle observables 

Linear responses cannot in general be expressed in 
terms of only two-particle observables but they can nev- 
ertheless be calculated with the formalism of the TPEM 
without introducing random vectors, uncontrolled ap- 
proximations or systematic errors, as will be shown. The 
discussion here will be focused on the optical conductivity 
but the procedure can be generalized to other observables 
as long as there is an operator that plays the role of the 
current. The key step is based on a recent development 
by Weiss o 18 ' 19 in the context of the Anderson model. 

For any configuration of classical fields 0, the opti- 
cal conductivity^ can be rewritten as (/(x) is the Fermi 
function): = ^ £ AA , | (A| J;c |A') | 2 (/(/3(£a' -fi))- 

f(/3(Ex-(i))5(uj-Ex+Ex>), where |A) and E x denote the 
eigenvectors and eigenvalues of the one-electron Hamil- 
tonian. Note that E\ is related to the normalized Hamil- 
tonian eigenvalues simply by E\ = ae\ + b (see Eq. (Q). 
Now, following Reference fl8L ct{uj) can be expressed in 
terms of the function j^(x,y): 

U(x,y) = ^^|<A|^|A0| 2 5(.T-e A )%-eA') (11) 

AA' 
1 f°° 

o>(w) = — j^(x,x + u))[f(P(ax + b- fji))- 

10 J-oo 

f(/3(ax + b- fi + uj))]dx. (12) 

The function j^(x, y) can be expanded in Chebyshcv 
polynomials: j<p{x,y) = YlmnVmni^TmixjT^x), with 
r, mn {cj>) = Tr(T m (H ((/>)) j x T n (H ((/>)) j x )/N. If the above 
trace is evaluated in the real-space basis, {\i >, er}, of the 
one-electron sector, then it is possible to express rj mn in 
terms of the coefficients fiia-ja'-.m defined in Eq. ||SJ|: 

Vmn(4>) = ^ ^2 ViviJ<?2\m((l))(jcr2\jx\k<T 3 ) X 

l,j,k,l <T\ ,<T2 j(T3 ,(74 

XPk<T3,la*in(.4>)(l< T 43x\ i<T l)> ( 13 ) 



E. Parallel Computation 

Another advantage of the TPEM, is that it can be par- 
allelized effectively up to 50-100 processors, because the 
calculation of the moments for each basis ket \v) is inde- 
pendent. Thus, the basis can be partitioned in such a way 
that each processor calculates the moments correspond- 
ing to a portion of the basis. It is important to remark 
that the data to be moved between different nodes are 
small compared to calculations in each node. Indeed, the 
communication time is proportional to MNpe and com- 
munication among nodes is mainly done here to add up 
all the moments. The exact diagonalization algorithm 
can be parallelized only up to very few processors even 
when considering sparse matrices. 

The calculation of observables that will be presented 
in the next section was carried out by using parallel com- 
putation. The complexity scales with the number of pro- 
cessors up to the size of the one-electron Hilbert space, 
Ndim- However, the observable that needs to be calcu- 
lated more frequently is the difference in effective action. 
Due to the truncation in this action difference that takes 
place when the TPEM is used (see e. g. Ref. 0), the sum 
in this case is performed only up to the size of the trun- 
cated Hilbert space, A. This implies that the inverse 
complexity scales with the number of processors, Npe 
only up to Npe = A and not Npe = Nd im . Although 
this imposes a limitation of the parallelization algorithm, 
since in general A -c Ndim, energy integrated observables 
usually converge for a cutoff M < 20 and less than A pro- 
cessors is enough in most cases 

However, energy dependent or many-body observables, 
e. g. the density-of-states and the optical conductivity, 
need more moments (usually 100 or more as will be dis- 
cussed in the next section) for their calculation within 
the same error. Now, increasing the number of moments 
M implies an increase of the truncated Hilbert space of 
size A. But even though now the number of moments has 
to be increased, the threshold for the scaling of the paral- 
lelization algorithm has also increased. Therefore, these 
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types of problems are particularly suited to be studied 
efficiently with supercomputers with thousands of nodes. 



III. APPLICATION: DOUBLE EXCHANGE 
MODEL FOR MANGANITES 



A. Microscopic Model 
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FIG. 1: (a) Kinetic energy of a single configuration of classi- 
cal spins on a 8x8 lattice for model Eq. 1151 calculated with 
the DIAG and TPEM algorithms. The configuration was ob- 
tained after evolving the model with Monte Carlo with 600 
steps for Jh = 4t and p = — At, T = Q.02t and Jaf = 0. 
t pr = 10 -5 and €tr = 10~ 6 . (b) Density-of-States for the 
model for the same configuration and parameters as in (a). 
Algorithms and different expansion cut-offs are indicated. 

Manganese oxides, known as manganites, are 
correlated-electron systems of much current interest in 
condensed matter physics. This interest started with the 
discovery of the so-called "colossal" mangetoresistancc 
(CMR) effect, an effect three orders of magnitude larger 
than the one observed previously on superlattice films^i 

Considerable progress has been achieved in under- 
standing many properties of these compounds by using 
double exchange models^ The electron-phonon coupling 
is usually neglected but it can and has been included; 
the phonon degrees of freedom are, then, an additional 
classical field. In the ferromagnetic phase the electrons 
directly jump from manganese to manganese spin and 
their kinetic energy is minimized when these spins are 
aligned. Moreover, the unbiased study of spin-fermion 
models for manganites as presented here can produce 
other long- and short-range ordered phases, e. g. antifer- 
romagnctism, inconmensaturate phases (to be described 
later) and phase separation. In general, the correct 
ground states and very good estimates of Curie and Neel 
temperatures have been obtained with the DIAG method 
applied to double-exchange models^. The Hamiltonian 
of the system can be written asS2i24: 



— Jh 22 Si-Si 



-J AF J2SrS 3 , (15) 



where cJ CT creates a carrier (e g electron) at site i with 
spin a. The carrier-spin operator interacting ferro- 
magnctically with the localized Mn-spin S 1 , is Bi = 
J2 a a c\ a &a,i3Cif3- For the manganites, the Hund coupling 
Jh is large compared to the hopping energy scale t and, 
as a consequence, suppresses double occupancy. In this 
respect, a large Jh coupling acts as a large Hubbard 
repulsion. The last term is the superexchange between 
localized t2 g spins. 

Static properties of this model, e. g. the local spin 
magnetization and energy of the system have been pre- 
viously studied with the TPEM on large lattices^ Here 
the focus will be on more complex observables. Before 
showing results for dynamical observables we will present 
first results for a two-particle static observable, namely, 
the kinetic energy, K. A simulation was performed in 
the ferromagnetic region of Hamiltonian Eq. (|15[) with 
Jh = 4i and quarter filling. The kinetic energy is shown 
in Fig. for a fixed configuration of classical fields 
(spins, Si, in this case) calculated both with the diago- 
nalization method and with the TPEM for various values 
of the cutoff, M. The configuration of classical spins was 
selected after evolving the system with Monte Carlo at 
T = 0.02i. K converges to the exact value within an 
error window of less than 2% for M > 100 which implies 
that the required cutoff is larger than for one-particle or 
classical observables ^ 

Periodic boundary conditions will be used throughout 
the following calculations. The unit of energy is the hop- 
ping parameter, t. 



B. Density of States and Spectral Function 
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FIG. 2: DOS at three temperatures on a 20x20 lattice cal- 
culated with the TPEM and the Monte Carlo algorithm, 
M = 200, for model Eq. at J H = St, (n) = 0.5 for (a) 
Jaf = and (b) Jaf — 0.02. The location of the chemical 
potential is indicated by the vertical dashed line. 
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FIG. 3: Spectral function for the double exchange model 
calculated on a 20x20 lattice with the TPEM at J H = St, 
n = 0.5, T — 0.02t (FM system). The curves are parameter- 
ized in terms of (k x ,k y ) from (0,0) to (7r,0) to (n,n). The lo- 
cation of the Fermi energy is indicated by the vertical dashed 
line. There is a finite DOS at the Fermi energy. The dashed 
fine is the spectrum calculated analytically for a perfect fer- 
romagnet. 




FIG. 4: For Jaf = 0.15 and the same parameters as in Fig.|3] 
corresponding to an insulator, A(k,iu) shows a clear gap and 
a more complex form. The vertical dashed line indicates the 
position of the Fermi energy. 



We will analyze Eq. (fPo|> at quarter filling ((n) = 0.5) 
and Jh = 8i. For Jaf = the ground state of the system 
is ferromagnetic. As Jaf is increased the superexchange 
coupling competes with the double-exchange term and 
there exists a critical Jaf beyond which the system's 
ground state is no longer ferromagnetic. This critical 
Jaf was found to be 0.05 ± 0.01 in two dimensions and 
0.12±0.01 in one dimension. Above this value, the struc- 
ture factor presents peaks at incommensurate momenta 
(and if Jaf is large enough the system becomes antifcr- 
romagnctic). In the next section we will show that the 
ferromagnetic region is metallic, whereas the region with 
incommensurate momenta is insulating. 



First we show that the DOS calculated with the TPEM 
gives reasonable results by comparing with the ones ob- 
tained with the DIAG (exact) method. This is presented 
in Fig. [iy> for the same parameters as in Fig. QJi- The 
finite size effects of the DOS are only reproduced by the 
TPEM with 100 moments or more. Unlike energy inte- 
grated quantities, energy dependent quantities need more 
moments for convergence, as expected. The DOS in the 
ferromagnetic phase has two bands centered at approxi- 
mately ± Jh each with bandwidth equal to 4t. 

Now that we know that the TPEM yields reliable re- 
sults for the DOS we proceed to perform Monte Carlo 
simulations. In Fig. |2 the DOS is shown for differ- 
ent temperatures on a 20x20 lattice at (n) = 0.5 and 
Jh = 8t. The dependence with temperature is very mild 
at Jaf = (Fig. and there is a slight decrease of 
the DOS at the Fermi energy for J A f = 0.02 (Fig. 
which is still ferromagnetic. There is no indication of 
a pscudogap forming at higher temperatures in the fer- 
romagnetic region. However, in photoemission studies 
for manganites^ the DOS of metallic samples was found 
to depend more strongly on temperature and to display 
a more complicated structure than the one shown in 
Fig. |21 It is very likely that the inclusion of phonons 
or many orbitals are necessary to be able to reproduce 
these features. The influence of these corrections will be 
prominent mainly near the critical Jaf that separates the 
metallic from the insulating regions. This suggests that 
the relevant regime of our model for typical ferromagnetic 
samples lies very near the metal-insulator transition that 
occurs with increasing Jaf- We will see evidence of this 
also when we examine the optical conductivity and con- 
ductances in the next section. 

Our tests show that the calculation of the spectral 
function with Eq. IjlOfl requires at least M = 100 for 
this model. Using a large value for M, A(k,u) is plotted 
in Fig.|3and^on a 20 x 20 lattice for Jh = 8i at low tem- 
perature (T = 0.02i) and quarter filling. A total of 1000 
thermalization steps and 1000 measurement steps were 
done first with M = 40. An additional 100 steps were 
used to calculate A(k, to) with a large number of moments 
for the expansion. For Jaf = 0, Fig- El corresponding to 
a ferromagnet, there is a single peak for every value of k, 
as expected since the system is ferromagnetic. There is 
also a clear Fermi surface satisfying cos(fc x ) + cos(fc. y ) = 
and the peaks follow the dispersion calculated analyti- 
cally for a perfect FM. But for Jaf = 0.15 in Fig. H 
A(k, ui) presents a complex form and also a gap opens at 
the Fermi energy. It is worth noting that the structure 
factor at this value of Jaf has peaks at incommensurate 
momenta, e. g. (k x , k y ) = (tt, 0) that are precursors of the 
antifcrromagnctic phase that forms at even larger Jaf- 

In all cases, the moments where multiplied by damp- 
ing factors to reduce the high frequency oscillations as 
explained before. 
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C. Optical Conductivity 

The optical conductivity in units of e 2 /h is shown in 
Fig. [3] on a 6 3 lattice for two configurations of classical 
spins obtained at two different temperatures. The figure 
shows that the TPEM with 200 moments gives results 
consistent with the diagonalization method. The value 
of Tq — 0.09 ± 0.01 was determined from the maximum 
distance correlation between classical spins (Fig. [HJ) taken 
as distance — Ly/3/2, with L — 6 being the linear length 
of the lattice. Below Tc there is a single peak at to = 
and above Tc there is an additional peak at u) = 2Jh that 
increases in intensity as temperature is increased. This 
peak at finite lu ~ 2Jh is due to transitions from both the 
spins up and down since they have similar occupations 
above Tc- On the other hand, below Tc the electronic 
states with spin antiparallcl to the t2 g spins have high 
energies and, therefore, the transition across the 23 h gap 
weakens and only the "Drude" peak at ui = is left. 

Fig-0is in qualitative agreement with the optical con- 
ductivity spectra of Lai-^Sr^MnOa (see e. g. Ref. l2r|) 
where the intensity is transfered from low to high en- 
ergy as temperature increases. Then, our three dimen- 
sional results suggest that the spin-fermion model con- 
sidered here contains the basic physics to reproduce the 
behavior of the optical conductivity with temperature at 
least in the metallic regime. A more accurate treatment 
would include many orbitals, in particular if orbital or- 
dered phases, e. g. the CE phased, are relevantj^&Si. 
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FIG. 5: Optical conductivity on a 6 lattice calculated with 
the TPEM and DIAG for the same configurations of classical 
fields at Jh = 8t and quarter filling, (a) T = 0.02t calculated 
with the DIAG (b) T = 0.02t calculated with the TPEM, (c) 
T — 0.2t calculated with DIAG and (d) T = 0.2i calculated 
with the TPEM. In (b) and (d) M = 200. 

Now the TPEM will be used to determine the metal 
to insulator transition that occurs in the model when 
a Jaf 7^ is considered in Eq. 1)15(1. This transi- 
tion was found in one dimensional finite chains at zero 
tcmperaturaSi. The quantity of interest is the limit of 
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FIG. 6: Maximum distance spin-spin correlation, C(d max ), 
vs. temperature for the ferromagnetic regime (Jaf — 0) on 
6 3 and 8 3 lattices. d max = \/3L/2 where L = 6 or L = 8. 



ct(uj) as u — ► 0, yielding the conductance C. A total 
of 1000 Monte Carlo steps for thcrmalization and 1000 
measurement steps were performed with AI = 40. An 
additional 100 steps were used to calculate dynamical 
properties with a large number of moments for the ex- 
pansion (usually M = 200). We will examine the three- 
dimensional case for two values of the superexchange cou- 
pling Jaf- In Fig. [7|i for Jaf — (corresponding to 
a ferromagnetic ground state), C is plotted as a func- 
tion of temperature on 6 3 and 8 3 lattices with Jh = 8t 
and [i — — 8t, i. e. for a density approximately equal 
to quarter filling. As expected C decreases monotoni- 
cally with increasing temperature. On the other hand, 
for Jaf = O.lt the ground state is not ferromagnetic 
anymore but it is given by states with incommensurate 
momenta instead. These states are insulating as can be 
seen from Fig. {7]p. Therefore, the present method allows 
for the analysis of the metal to insulator transition with 
increasing Jaf in the three dimensional model. It was 
essential to be able to run simulations on 6 3 and 8 3 lat- 
tices since 4 3 is known to show strong finite size effects 
(e. g. the low temperature magnetization increases with 
temperature in that case). Simulations on 8 3 or more 
lattice sites are not possible with the DIAG method. 

Note that the standard dependence of the resistivity on 
temperature for CMR manganites cannot be reproduced 
with the model given by Eq. I|15|l : the calculated resistiv- 
ity has no peak at finite temperature, only a monotonous 
behavior. However, there is indication^ that chemical 
disorder could transform an insulator into a poor metal 
at low temperatures. For this to happen the parameters 
of the system (Jaf in our model) have to be located very 
near the border of the metal-insulator separation. This 
is in agreement with the experimental observation^ that 
CMR manganites in the ferromagnetic region have poor 
conductance and with the results presented before for 
the density-of-states and spectral function and provide 
further evidence that this is the relevant regime to in- 
vestigate for possible explanations of the CMR effect. In 
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other words, if small corrections, e. g. disorder, many or- 
bitals or weakly-coupled phonons, are included far from 
the critical Jaf, they will not affect the qualitative trends 
shown here but they are expected to be very important 
near this critical Jaf- 

The Hamiltonian parameter corresponding to the su- 
perexchange interaction between t2 g spins, Jaf, drives 
the three-dimensional system from a metal to an insula- 
tor. In the phenomenological model of Eq. (|15(l . Jaf is 
an effective interaction that competes with the double- 
exchange mechanism. This competition is related to 
the tolerance factor of the samples, or the deviation of 
the Mn-O-Mn angle. Therefore, calculations at differ- 
ent values of Jaf can simulate the transport and mag- 
netic properties of different compounds. For instance, 
in the Pri_ 2; (Cai_j,Srj,) a ,Mn03 system^ distortions in 
the MnC>6 octahedra weaken the double-exchange mech- 
anism and the superexchange interaction becomes impor- 
tant as well as collective Jahn- Teller distortions. Even 
though in this work we have included only the superex- 
change coupling, the qualitative consequences of this 
competition can be seen in Fig. Moreover, with the 
formalism presented in this work for the TPEM it will be 
possible to obtain linear responses of Hamiltonians that 
include Jahn- Teller phonons and many orbitals which is 
not the case with the DIAG method because the Hilbcrt 
space of these systems is at least twice as large as the one 
for the one-orbital model, i. e. the computations would 
be 16 times slower. The complexity of the TPEM scales 
linearly with the size of the one-electron sector. 
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FIG. 7: cr(0) vs. temperature at Jh 
on 6 3 and 8 3 lattices for (a) Jaf = 
showing metallic and insulating behavior with temperature 
respectively. The results were obtained with the TPEM and 
M = 200. 



IV. CONCLUSIONS 



The TPEM was extended to include the calculation 
of two- and four-particle observables without introduc- 
ing systematic errors or uncontrollable approximations. 
By comparing with the exact algorithm it was found that 
energy-dependent observables need more moments than 
energy integrated observables but the results converge 
with M « 100 moments. The formulas were applied 
to a spin-fermion model for manganites on a 20x20 lat- 
tice to analyze the spectral function for the ferromag- 
netic and incommensurate phases. We concluded that 
the simple model of Eq. I|15|l cannot reproduce the DOS 
of typical CMR samples. The inclusion of chemical disor- 
den or weakly-coupled phonons would be necessary and 
these corrections should be important near the border 
that separates the ferromagnetic from the region with 
incommensurate momenta in the phase diagram. Similar 
conclusions were found by studying the metal to insu- 
lator transition with increasing superexchange coupling 
that occurs in the three dimensional model. 

The truncated polynomial expansion allowed us to cal- 
culate dynamical properties of a spin-fermion model for 
manganites on up to 8 3 lattices without approximations. 
Even larger lattices will be reachable with the new gen- 
eration of supercomputers that can accommodate easily 
thousands of nodes, as explained before in EIII El Big 
lattices are particularly necessary in simulations includ- 
ing chemical disorder. Moreover, systems with realistic 
band structures^ will require an increase in the size of 
the Hilbert space and the calculation of dynamical ob- 
servables in these more realistic models will benefit from 
the present methods. 
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